Full counting statistics for transport through a molecular quantum dot magnet 
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Full counting statistics (FCS) for the transport through a molecular quantum dot magnet is 
studied theoretically in the incoherent tunneling regime. We consider a model describing a single- 
level quantum dot, magnetically coupled to an additional local spin, the latter representing the 
total molecular spin s. We also assume that the system is in the strong Coulomb blockade regime, 
i.e., double occupancy on the dot is forbidden. The master equation approach to FCS introduced 
in Ref. [12 is applied to derive a generating function yielding the FCS of charge and current. In 
the master equation approach, Clebsch-Gordan coefficients appear in the transition probabilities, 
whereas the derivation of generating function reduces to solving the eigenvalue problem of a modified 
master equation with counting fields. To be more specific, one needs only the eigenstate which 
collapses smoothly to the zero-eigenvalue stationary state in the limit of vanishing counting fields. 
We discovered that in our problem with arbitrary spin s, some quartic relations among Clebsch- 
Gordan coefficients allow us to identify the desired eigenspace without solving the whole problem. 
Thus we find analytically the FCS generating function in the following two cases: i) both spin 
sectors (J = s ± 1/2) lying in the bias window, ii) only one of such spin sectors lying in the bias 
window. Based on the obtained analytic expressions, we also developed a numerical analysis in 
order to perform a similar contour-plot of the joint charge-current distribution function, which have 
recently been introduced in Ref. Il3t here in the case of molecular quantum dot magnet problem. 



I. INTRODUCTION 



Molecular electronics has emerged as one of the promis- 
ing subfields of nanophysics. Individual nano objects 
such as molecules and carbon nanotubes can be con- 
nected to leads of different nature, and the degrees of 
freedom of molecules could be used in principle to mod- 
ulate quantum transport between the leads [l], 0, H| : it is 
hoped that these degrees of freedom could allow specific 
functions of nanoelectronics. A large body of work has 
focused on the potential applications of the specific en- 
ergy spectrum of molecules, as well as their vibrational 
degrees of freedom. [1] Recently, there has been some 
focus on transport through molecules which have their 
own spin, for instance because a specific atom of this 
molecular system possesses a magnetic moment. This is 
unlike "usual" molecular quantum dots where the molec- 
ular spin corresponds to the effective spin of the itinerant 
electron, giving rise for example to Kondo physics. [5[ 
There are many experimental proposals for such molec- 
ular magnets. Mni2 acetates or carbon fullerenes with 
one or more atoms trapped inside are examples. Depend- 
ing on the sample, and on the structure of the molecule, 
the spin orientation can have a preferred axis or it can 
have a more isotropic nature. There have been recent ex- 
periments on similar systems, with applications to spin 
valves. [f| 0] In this paper we focus on the isotropic 
regime. As in Ref. 8], where the current through an 
cndohedral nitrogen doped fullcrcnc was computed, we 
describe transport in the incoherent tunneling regime, 
and proceed to compute the full counting statistics for 
transport through a molecular quantum dot where the 
dot electrons have an exchange coupling with an impu- 
rity spin. 



Indeed, in mesoscopic physics, [1(3] the current volt- 
age characteristics alone is in general not sufficient to 
fully characterize transport. The current fluctuations 
can provide additional information on the correlations 
in molecular quantum dot, but so do the higher current 
moments, or cumulants, such as the "skewness" (third 
order) and the "sharpness" (fourth order). Such higher 
order cumulants add up to the so-called full counting 
statistics (FCS). [ll[ Yet here we are interested in ob- 
taining the FCS of such devices using the master equa- 
tion approach. [12| We will also take the opportunity to 
present in detail the general framework for studying joint 
distributions of electron occupation of the dot and of cur- 
rent. A recent work by one of the authors [H[ applied 
to a single level quantum dot shows that there are clear 
correlations between current and occupancy of the dot: 
there, the quantum fluctuations caused by increasing the 
tunnel coupling to the leads has a definite impact on the 
number distribution. Here we will extend these concepts 
to the case of a molecular quantum dot with an impurity 
spin. 

The FCS has remained until recently a highly theo- 
retical concept, with few experimental applications. The 
third moment has been measured for a tunnel junction, 
[l4l , [l5j as the effects of the measuring apparatus have 
been pointed out. [la ] For the regime of incoherent trans- 
port through quantum dots two outstanding experiments 
have been performed: [13, a point contact is placed in 
the close vicinity of the quantum dot, and by monitoring 
the current in this point contact, it is possible to obtain 
information on the time series of the current in the dot. 
From this time series, higher moments of the current can 
be computed. These recent experiments thus provide ad- 
ditional motivation to study the FCS in this novel system 
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of a molecular quantum dot with an external spin. 

This paper thus aims at calculating the FCS of trans- 
port through a molecular quantum dot magnet, under 
the following circumstances: 

1. The quantum dot is in the strong Coulomb block- 
ade regime, i.e., the system is subject to strong 
electronic correlation. 

2. The coupling between the dot and the reservoirs is 
relatively weak, and tunneling of electrons between 
them can be considered incoherent (incoherent tun- 
neling regime). 



merical analysis to perform a contour-plot of the joint 
charge-current distribution function. 

The paper is organized as follows. The first two sec- 
tions after Introduction are devoted to the preparatory 
stage, i.e., Section II describes our model for the isotropic 
molecular quantum dot magnet, while Section III reviews 
the master equation approach to FCS. Then, in Section 
IV the analytic expressions for the FCS generating func- 
tion will be obtained, with the help of quartic identities 
among Clebsch-Gordan coefficients. Section V describes 
the contour plot of joint distribution functions, before 
the paper is concluded in Section VI. 



3. The molecular spin is isotropic, and coupled to the 
spin of conduction electron through a simple ex- 
change interaction. 

Assumption (2) allows us to employ the master equation 
approach for studying the transport through molecular 
quantum dot. As assumption (1) forbids double occu- 
pancy on the dot, we consider the cases where the dot is 
either empty, or occupied by only one conduction elec- 
tron. 

When the dot is empty, the molecular spin s, together 
with its z-component s z , is naturally a good quantum 
number, whereas if the dot is occupied by one conduction 
electron, then it becomes the total angular momentum 
j = s + a which is a good quantum number on the dot. 
Let us now follow the time evolution of the dot state, 
e.g., using a master equation. The number n of con- 
duction electrons on the dot fluctuates in time between 
n = and n = 1. Such transitions are controlled by the 
so-called rate matrix L, which will be more carefully de- 
fined in Sec. II. This rate matrix L can be decomposed 
into 2x2 = 4 charge blocks , corresponding to two charge 
sectors n — 0,1. Then each block is further character- 
ized by either the molecular spin s z or the total angular 
momentum j (and j z ) depending on whether n = or 
n = 1. Naturally, off-diagonal blocks of the rate matrix 
describe transition between different angular momentum 
states. The rate of such transitions is proportional to 
the square of a Clebsch-Gordan coefficient between the 
initial and the final angular momentum states. 

In the master equation approach to FCS, the deriva- 
tion of FCS generating function reduces to solving the 
eigenvalue problem of a modified master equation with 
counting fields. By using an algebraic identity, this eigen- 
value problem can be projected down to a slightly more 
complicated problem within the n = charge sector, 
which can be fully characterized by the molecular spin 
s z . Our main discovery is that some identities among 
Clebsch-Gordan coefficients allow us to identify the de- 
sired eigenspace without solving the whole problem. In 
this reduced problem Clebsch-Gordan coefficients appear 
at quartic order, since off-diagonal blocks are projected, 
together with the n = 1 charge sector, down to the n = 
subspace. Based on the analytic expressions for the FCS 
generating function thus obtained, we also developed nu- 



ll. MOLECULAR QUANTUM DOT MAGNET 

We consider a molecule connected to two electrodes, 
which are specified by an index \ = ±1. \ — 1 corre- 
sponds, e.g., to the left (x = —1 to the right) electrode. 

We focus on the strong Coulomb blockade regime, as- 
suming that the number of electrons which can be added 
to the dot is restricted to one. Because the extra-charge 
electron state on the dot carries a spin, it is locally cou- 
pled to the impurity spin on the molecule via an exchange 
coupling Hamiltonian. 



A. Model Hamiltonian 

The total Hamiltonian H of the system consists of two 
parts: (i) the Hamiltonian Hq of the isolated parts, and 
(ii) the tunneling Hamiltonian H tU n, i.e., H = H + H tun . 
H can be decomposed further into the dot and lead 
parts: H = X) x =±i -^liad + ^dot- For normal metal 
leads, the Hamiltonian of left (\ = 1) and right (x = — 1) 



leads reads, H, 



(x) 



Ejt-ei^^cl* 5 , where c[ x J aum- 



Mcad L^ko c fe ^ka u fc<7 ' ^ka 

hilates an electron with momentum k, spin a and en- 
ergy in the lead \- The tunnel Hamiltonian reads: 



H 



tun = Y<kax T kcr C ka d ° + n - c -> wnerc tnc operator d 
(d)) annihilates (creates) an electron on the single dot 
level. The dot Hamiltonian reads: 

#dot = £dot 5Z n<T + Un l n l + #spin, (1) 
cr 

where n a — d\d a is the number of electrons on the dot 
level with spin a, and U is a charging energy contribution 
related to the capacitance with respect to the leads and 
to the gate. It will be assumed to be large in the present 
work, in order to prohibit double occupancy. Finally, 
H sp in is the coupling of the dot electron spin a with the 
impurity spin s. This is finite, by definition, only when 
there is an electron on the dot: 



H s , 



(n = 0) 

ga ■ s (n = 1) ' 



(2) 



where g is the exchange coupling constant. As expected 
the eigenstates of this contribution to the Hamiltonian 
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are the total angular momentum states of j = a + s. This 
is seen by rewriting the exchange term as a ■ s = {j 2 — 
s 2 — <t 2 )/2. The n = sector is characterized simply by 
the spin s of the magnetic impurity, and has 2s + 1 degen- 
erate states, \n — 0, s,s z ). When n = 1, the coupling g 
between the impurity and the dot electron spins requires 
to diagonalize the interaction with the total angular mo- 
mentum states, i.e., \n = l,j = s±l/2, j z ). The two spin 
subsectors, j = s + 1/2 and j = s — 1/2 correspond to an 
energy eigenvalue of, respectively, e s +i/2 = £dot + 5^/2 
and e s _!/ 2 = £dot — g(s + l)/2. The two subsectors each 
have a degeneracy of 2s + 2 and 2s, respectively. 

In our model Hamiltonian, the charging effects are 
taken into account via the bias- voltage dependence of the 
position of the dot level e^ot with respect to the chemical 
potentials of the leads, hl,r- A gate voltage V g and gate 
capacitance C g can be included in order to move the dot 
levels up and down. Here we focus on the case where the 
bias voltage always dominates on the temperature, so for 
later purpose it is sufficient to specify which level and 
how many of these are included in the bias window. 



B. Incoherent tunneling regime - the master 
equation 

The electron dynamics of transport through a molecu- 
lar quantum dot in the incoherent tunneling regime can 
be described by a master equation. [20( In such sys- 
tems as our molecular quantum dot Hamiltonian, the 
master equation describes a stochastic evolution in the 
space spanned by the eigenstates of the non-perturbed 
Hamiltonian Hq. Such states are, of course, of quan- 
tum mechanical nature, but the evolution of the system 
is governed by a deterministic stochastic process. Jumps 
between different states are to be stochastic events, and 
are also to be Markovian (no correlation between suc- 
cessive tunneling events). Such a situation is justified in 
the so-called incoherent tunneling regime, where the the 
non-perturbed Hamiltonian Hq is only weakly perturbed 
by the tunneling Hamiltonian iftun , an d the tunneling of 
electrons through the dot is still considered to be sequen- 
tial. 

A simplification specific to our model is that the lead 
degrees of freedom can be traced out, and the master 
equation describes actually the time evolution of the dot 
states. We will often denote a (quantum mechanical) 
dot state \n,j,j z ), simply as a in the master equation. 
The probability p a {t) with which the system, i.e., our 
molecular quantum dot, finds itself in a state a at time 
t satisfies, 

^=7«Pa(*) + £r a/jP/3 (t). (3) 

p 

T a p is the transition probability for a jump, a <— /3, 
which can be calculated, e.g., from Fermi's golden rule. 
The subscripts a and f3 take all the possible values, 



ax, ct2, • ■ • in the space of physically relevant states. Con- 
servation of probability imposes that the sum of all the 
components in a given column (5 of matrix L vanishes. 
Thus 7 Q = — J2p T/3q denotes the decay rate of a state 
a. In matrix notation, this reads dp(t)/dt = Lp{t), where 
L = 7 + r, is called the rate matrix. 7 and T are, re- 
spectively, the diagonal and off-diagonal parts of the rate 
matrix L. 

In order to apply Fermi's golden rule for evaluating 
T a p, one must go back to the Hamiltonian, Hq of the 
isolated parts, and one consider the transition rate Tfi 
associated with a tunneling perturbation H tun , from an 
initial state \i) to a final state |/): 

2-7T 9 

r fi = T \(i\H tun \f)\ 2 s{E i -E f ) .. (4) 

The eigenstates \i),\f) are many electron states, but 
they factorize into a direct product of many-electron lead 
states and a dot state. Therefore, after taking into ac- 
count all the lead degrees of freedom, the above transition 
rates are rewritten in the form of transition rates between 
different dot states a, 0, i.e., T a/ 3 for a transition, a <— (3. 

We also make the standard assumption that the hop- 
ping amplitude T^' depends weakly on energy for states 
close to the Fermi levels of the leads. For normal 
leads it is also spin independent. Thus Fermi's golden 
rule gives basically a bare tunneling rates, defined as 
Fl,r = 2nvTl R , weighted by a Clebsch-Gordan coef- 
ficient squared, which we will discuss later, v is the (con- 
stant) density of states at the Fermi level in the leads. 
We also use the explicit notation L, R rather than \ f° r 
specifying the leads. 

Before giving explicit matrix elements of the rate ma- 
trix L, it is convenient to decompose it into four blocks, 
corresponding to the two charge sectors n = 0, 1 as 
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The off-diagonal block matrices B and C are rectan- 
gular matrices of size, (i) (2s + 1) x (4s + 2) and (ii) 
(4s + 2) x (2s + 1), respectively. They describe tunneling 
of an electron either (i) out of the dot into one of the 
reservoirs, or (ii) onto the dot from one of the reservoirs. 
Both have two spin sectors, corresponding to two possi- 
ble values of the total angular momentum j = s ± 1/2, 
i.e., B = {B s _ l / 2 ,B s+1 / 2 ) and C = (C s _ 1 / 2 , C S+1 / 2 ) T , 
where the C T is, for example, a transposed matrix of C. 
The (s z , jz)-component of the submatrices B s ±x/2 and 
the {j z , s 2 )-component of C s ±x/2 are given, respectively, 
as 

Bj{s z ,j z ) = 

[T L {1 - f{ej - ti L )} + T R {1 - f(ej - m)}] 
x\{s,s z ;l/2,a z \j,j z }\ 2 (6) 
Cj{jz,S z ) = 

[ri/(ej - (i L ) + T R f{ej - /j, R )] 
x\{s,s z ;l/2,a z \j,j z )\ 2 , (7) 
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where j = s± 1/2, j z = -j, ■ ■ ■ ,j. 

The diagonal blocks A and D are themselves diagonal 
matrices, because the dot Hamiltonian is diagonal in the 
chosen basis. Their elements are specified by recalling 
that the sum of all the components in a given column of 
L vanishes identically, as a result of the conservation of 
probability. 

The matrix D has two spin sectors, D = 
diag(Z3 s _ 1 / 2 , D s +i/2)- The two diagonal blocks Dj are 
also diagonal matrices, Dj(j z ,j' z ) = Dj(j z , j z )S(j z , j' z ). 
Let us focus on the column of L belonging to the sector j 
and j z . The conservation of probability associated with 
this column requires, 

D 3 (jz,jz) = ~ B j( s zJz) 

S z — — S , ■ ■ ■ . S 

= -T L [1 - f(e 3 - ml)] - r K [l - fa ~ » a )] 



The quantity of our eventual interest is the joint statisti- 
cal distribution function P(qi, q2, ■ ■ ■ ) of such observables 
averaged during time r, i.e., 

P(qi, <72, • • • ) - (<%i - Qi)<%2 - Q2) ■ ■ ■ >, (11) 

where (■ • • ) represents a stochastic average, whose mean- 
ing will be more rigorously defined in Sec. Ill B. The 
joint distribution function P{q\, q2, ■ ■ ■ ) is related to the 
generating function via a Fourier transformation: 



(8) 



i.e., Dj(j z ,j z ) actually does not depend on j z : this sim- 
plification is due to the identity (|A3j) . Similar quadratic 
relations among Clebsch-Gordan coefficients are also 
summarized in the Appendix. 

The diagonal elements of matrix A can be deter- 
mined in the same manner. Recall that A(s Zl s' z ) — 
A(s z , s z )S(s z , s' z ). This time one fixes s z in order to spec- 
ify a particular column. Then, in order to apply the con- 
servation law, we sum over the two spin sectors j and 
their associated j z , to find, 



A(s z ,s z ) = - Y E Cj(j z ,s z ) 

j=a±l/2jx=-j,— ,3 

2s r 

r i/(e s -i/2 - Ml) + r fl /(e s _ 1/2 - fj, R ) 



2s + 1 
2s + 2 r 



2s + 1 



r L/(e s+ i/ 2 - Ml) + IW(e s+1/2 - fi R ) .(9) 



Here we used the identity (|A1[) in the Appendix. 

The explicit matrix elements of the rate matrix are 
thus given. The discussion now turns to the description 
of FCS in the framework of master equation. 



III. MASTER EQUATION APPROACH TO FCS 

We are interested in the statistical correlations of such 
observables as charge N(t) or current I(t), measured dur- 
ing a time interval to < t < io+T. In the master equation 
approach, an observable Qf(t) (/ = 1,2, ••• ;Qi(t) — 
N(t),Q2(t) = I(t),---), is introduced as a stochastic 
variable which evolves as a function of time through the 
master equation. On the other hand, what is experimen- 
tally relevant is their time average during the measure- 
ment time r, 



yto+r 

Qf = ZI dtQ f {t). 



1 



E 

9l>92," 



P(qi,q2,---)exp Y 

/=1,2,- 

= exp [fi($i,$ 2 , •••)]• 



where f2($i, $2, ■ • • ) is called the cumulant generat- 
ing function (CGF). The (ni, 712, ■ • • )-th order cumu- 
lant C ni:n2l ... is deduced from the generating function, 
fi($i,$ 2 ,---) as, 



Cm ,ri2 



(-iy 



$ 2 , 



*1,*2, 



(10) 



Note that the CGF, fi($i, $ 2 , ■ ■ • ) and the full distribu- 
tion function P(qi, q%, • ■ ■ ) contain the same information. 
In the following, we will focus on how to calculate such 
CGF in the framework of the master equation approach. 

The purpose of this section is to review the master 
equation approach to FCS developed in Ref. [IH, and 
adapt it to a more specific case of our interest. A close 
comparison with the time-dependent perturbation theory 
of quantum mechanics reduces the calculation of CGF to 
solving a modified master equation with counting fields. 



A. Two types of observables and their counting 
statistics 

In Sees. IV and V we discuss counting statistics of 
charge N and current /. Other observables, such as the 
total angular momentum J on the dot, its z-component 
J z , or the spin current I s through the dot, can be han- 
dled in a similar way, which we outline below. Note that 
we can make a distinction between two category of ob- 
servables: "charge-like" and "current-like" observables. 
The charge N (or J, ■ ■ ■ ) is a property of the state of 
the dot, whereas the charge current, I c (or I 8 , ■ ■ ■ ) is as- 
sociated with transitions between different occupational 
states. The former has eigenvalues of a quantum me- 
chanical operator that can be simultaneously diagonal- 
ized with Hq, whereas the latter is related to the nature of 
the transitions. We may also classify the counting fields 
$1, $2, • • • into two categories: £ii £2, • • • and 771, 772, ■ ■ • ■ 

In order to give an unambiguous meaning to the 
stochastic average (• • • ) in Eq. (TT2"]) or (jTTJ) in the context 
of master equation, one has to go one step backward in 
regard to Eq. (fTU)) . i.e., one has to follow the time evo- 
lution of observables Qf(t). One may rewrite Eqs. (jTTj) 



5 



and jni as 



exp 
exp 



[n(ex,f2,- 



/ /' J 



(12) 



where the time integrals should be performed, e.g., from 
to to to + t. The stochastic average can be performed by 
giving more concrete expressions to the observables Nf (t) 
and If(t) using the language of stochastic process. 

Let us focus on a particular realization of the stochastic 
process by denoting {tj} the sequence of times at which 
jumps between different states occur, {aj} denotes the 
sequence of states such that the system stays in state aj 
during the period of tj < t < tj+i. A state \a) is specified 
by the quantum numbers n\,n-z, (e.g., ni = n, the 
charge, ni = j, the total angular momentum, n% = j z , 
its z-component, etc.). On the other way around, one 
denotes by ni(a), 712(0;), • • • the quantum numbers for 
specifying the state a. With this notation, Nf(i) can be 
written explicitly for a Markovian sequence ({tj}, {aj}): 



N f (t) = Y,n f (aj)6(t-tj)d(tj 



t), 



where 0(t) is the Heaviside function. 

A little more care is needed to give a similar expres- 
sion to current-like observables, reflecting the fact that 
they have a direction. In order to clarify this point, let 
us consider the case of our molecular quantum dot. The 
transition of the dot state is caused by tunneling of an 
electron onto or out of the dot. The change of the dot 
state before and after the transition is completely speci- 
fied by the change of the number of electrons An, and of 
the total spin Aj together with its z-component Aj z on 
the dot. However, in order to define a current associated 
with the transition, we need additionally an information 
on the direction of flow, i.e., an information on where the 
electron goes, to which lead, or from which lead. To mea- 
sure the current, one has to decide also where one mea- 
sures it. Such ambiguities are usually removed by con- 
sidering a net current which flows out of a specific lead, 
say, Xi i- e -i by considering A current li x ^ = An/ 

is induced by a tunneling, leading to the change of the 
dot state An/, such as An, Aj and Aj z , which occurs 
also through the lead %. This suggest in our Markovian 
sets ({tj}, {aj}) to add {xj} specifying through which 
lead the jump at time tj occurs. Once such Markovian 
sets ({tj}, {aj}, {xj}) are given, one can provide explicit 
formula to current-like observables: 



lf\t) = Y^^n f (t ] )6( X 



Xj)S(t-t j ), (14) 



where An/(i,-) = n/(oij+i) 
of transition at time t = ts. 



if(aj) specifies the nature 



Let us consider the case of our model again. We sep- 
arated the rate matrix L into two parts: the diagonal 
part 7 and the off-diagonal part T. We also introduced 
a decomposition into four charge sectors, which can be 
naturally applied to 7 and T as well: 



7 : 



j 7o 


° 


( Loo 





I 


71 J 


V 


Ln 









L01 


rio 


, 


1 = ( L 10 






r = 



Let us focus on the off-diagonal block matrices in T, T01 
and Tio- These two blocks correspond simply to dif- 
ferent An, i.e., Tio to An = 1, and Tqi to An = — 1. 
Such submatrices of T still contain contribution from dif- 
ferent leads, e.g., one can further decompose Tio into 
1 I^o . The whole T can be also decomposed 



with respect to x, i.e., T = X^T^, whereas one can 
divide the matrix T into smaller subsectors in terms of 
different An/. Thus very generally one can decompose 
the transition matrix T as 



r = 



E 

X,{An 



r 



(13) where {An/} = (An, Aj, • 



(x) 

{Any}' 



(15) 



B. Master equation with counting fields — 
prescription for obtaining the FCS 

Substituting Eqs. $T5$i and ^ into Eq. $FZ\) one finds 
the expression to be evaluated: 

exp Jo(£i,£ 2 ,--' ;vi,m,---) = 

n 

exp [iJ2^f ^2(t J+ i - tj)n f (aj) + 



/ 3=0 

n 

f 3=1 



(16) 



We need a prescription to compute the stochastic aver- 
age (■••)■ This is achieved by performing a perturba- 
tive expansion of the master equation with respect to 
off-diagonal components T. The resulting perturbation 
series is analogous to that of the time-dependent per- 
turbation theory of quantum mechanics. Note that the 
diagonal part 7 and the off-diagonal part T of the coeffi- 
cient matrix L = 7 + T generally do not commute. 

Because our observables are series of events which oc- 
cur successively on the time axis, the master equation is 
rewritten in a way which is analogous to the "interaction 
representation" of quantum mechanics, 



d , 
dt P ® 
Pit) 
f(t) 



f(t)p(t), 

exp [j(t - t )]p(t ), 

e -7(*-*o)p e 7(t-*o)^ 



(17) 
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Developing a perturbative expansion with respect to T, 
one finds the analog of the Dyson series for p(t) : 

P(t)=p(t )+ f dt't (t)p(f) 

OC p 

= X/ / dt 1 dt 2 ■ ■ -dt n 

n=0 Jt <t 1 <t 2 <---<t n <t 

xf(t n )---f(t 2 )T(t 1 )p(t Q ), (18) 

where we used p(to) — p(to)- Let us focus on the n-th 
order term in p(t) . Because Eq. (|18[) is written in matrix 
notation, we can further develop it as: 

f(t n )---f(t 2 )f(t 1 ) 

C*n — 1 »"' > Q 2 jCKl 

One can identify such terms as the realization of the 
stochastic process characterized by the sets {tj} and 
{a,}. Indeed, each of the n-th order term in Eq. fTS")) 
corresponds to a sequence, ao — > ai — > a2 — > ■ • ■ , with 
n jumps located at times {tj}. As the matrix 7 is the 
diagonal part of L, the system remains in the same state 
aj with probability e 7o, o(*i+ 1-t i) between the two jumps 
tj < t < tj + i, whereas the transition <x/-i — > aj occurs 
at time tj . Thus the probability for a sequence of events 
to occur is found to be, 



P(ao — > ai — > a2 

. . . p7a 2 (*3-«2)p 



->...) = 

The above identification in the framework of Eqs. (1 1 7[) 
and (|18p allows us to give an unambiguous meaning to 
the stochastic average of observables in Eq. (fT6|) . A 
prescription for calculating the average is given below. 

The first term in the exponential of Eq. p^|) consists of 
charge-like observables characterizing the state of the sys- 
tem for tj <t < tj+i. This contribution can be absorbed 
in the diagonal components of L, i.e., by the substitution: 

la ->1a(€l, &,•'•) = 7a + X ( 20 ) 

/ 

If one considers the distribution of charge in our molecu- 
lar quantum dot, this reduces simply to the replacement, 
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7o 
71 



l(0 = 



7o 
71 +£ 



The second term in the exponential of Eq. (fT6|) is, 
on the other hand, composed of current-like observables 
which are associated with transitions T a p between differ- 
ent states. In the perturbative expansion with respect to 
such jumps, these appear only at specific times t = tj. 
Eq. p9[) suggests that these contribution can be ab- 
sorbed, in the off-diagonal part T as follows: 



{An/} 



,(x) 

{An f } 



= , exp 

{An/} f 



/ 

(21) 



Note that only a part of T associated with the lead x 
which is involved in the measurement, i.e., is sub- 
jective to the shift rule. In our molecular quantum dot, 
r has the following simple structure, 







r (L) 
1 01 









r (R) 
1 01 





One also adopts the convention to measure the net cur- 
rent which flows out of the left reservoir onto the quan- 
tum dot. Then, the above replacement leads to, 











r {R) 
u 1 01 







The above arguments show that calculation of the CGF 
is actually equivalent to solving a stochastic problem, 
specified by the following modified master equation, 

j t P^mt) = £(£,i7)p(£, »?;*), (22) 
Ht,v) = 7(6,6, ••■) + r(T h , %,•••), (23) 

where p(£, n; t) is a vector: 

PaALmt) 

Its (£, ?y)-dependence is only written symbolically, 
but its implication would be clear: (£, 77) = 
(£i> 6> ' ' ' > Vii V21 " " " )• F° r a given initial condition 
Pi^iV'tto) = p(to), the CGF is simply related to p(£, r);t) 
as, 



exp 



The modified master equation (|23|) can also be integrated 
formally to give, 

p(£, m t) = exp (t - i )£(£, n) p(t ). 

Thus for a long measurement time r, the maximal eigen- 
value Am(C, 7 7) of the modified rate matrix L(^,n) domi- 
nates, i.e., 



(24) 



In the limit of vanishing counting fields: £, n — * 0, one 
can verify that Am(£, v) ~ * 0, the latter corresponding to 
the stationary state solution. 

We have thus established a recipe for obtaining FCS. 
We start with a master equation, Eq. ([3]), controlled by a 
rate matrix with its explicit elements calculated through 
the Fermi's golden rule. We then make the replacements 
(|20|21|) to define a new stochastic problem in terms of a 
modified master equation with counting fields: Eq. (|23p . 
Finally, our task reduces to finding the maximal eigen- 
value, Am(6 7 ?) °f the modified rate matrix L(^,rj) . The 
eigenspace corresponding to \m(£,v) is smoothly con- 
nected to the stationary state in the limit of vanishing 
counting fields. 
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IV. THE ANALYTIC SOLUTION 

The master equation describing an isotropic molecular 
quantum dot magnet in the incoherent tunneling regime 
was derived in Sec. II. A prescription for calculating FCS 
of transport through a quantum dot starting from such a 
master equation was given in Sec. III. The purpose of this 
section is to apply the formalism developed in Sec. Ill to 
the specific case of our molecular quantum dot magnet 
and obtain an analytic expression for the FCS of charge 
and current. 

In practice we have to solve an eigenvalue problem for 
the modified master equation with counting fields, which 
can be constructed by making suitable replacements of 
the diagonal and off-diagonal components as Eqs. (120]) 
and (|2ip . In this work, we will concentrate on the distri- 
bution of charge and current in our molecular quantum 
dot magnet. Our modified rate matrix is, therefore, a 
function of only two counting fields, say, £ and 77, i.e., 

L(£,r 1 )=\ (l) 70 ( R ) 01 + 01 (25) 

Its eigenvalues and the corresponding eigenvectors are 
also functions of £ and r\: 

LfovMS,v) = W,vWt,ri)' (26) 

The FCS of our problem is obtained as a cumulant gen- 
erating function of the charge-current joint correlation 
functions, which is identified as the maximal eigenvalue 
A (£,??) of the modified rate matrix, Eq. (|25|) . In the 
limit of vanishing counting fields ^, 77 — > 0, the eigenstate 
u(£,ri) corresponding to A(£, rj), collapses smoothly on 
the stationary state uq, i.e., 

lim rj) = u(0, 0) = uq, 
lim \(£,rj) = A(0,0) = 0. (27) 

In the same limit, Eq. (|2l)|) reduces to Eq. 

Obtaining FCS of a system thus does not require to 
solve completely the modified master equation, but sim- 
ply to identify the maximal eigenvalue, corresponding 
to the stationary state solution in the limit of vanish- 
ing counting fields. It is still rather surprising that one 
has an access to an analytic solution of such a problem 
for arbitrary s, because: 

1. The size of the matrix one has to diagonalize is 
"large", e.g., 6 x 6 for s = 1/2, and 3(2s + 1) x 
3(2s + 1) for a general s. 

2. For arbitrary s, one has to diagonalize such a ma- 
trix of that size, and it is generally not intuitive 
whether the obtained eigenvalues have a closed 
form as a function of s. 



For s = 1/2 the 6x6 eigenvalue problem, Eqs. (|26| and 
([23]) . can be treated analytically by reducing the size of 
the problem down to 2 x 2, using the identity, 



det \r n = det ( A ~ BD~ l C), (28) 



which is valid when det D/0. Here the four block matri- 
ces correspond to different charge sectors. One can apply 
the same identity to the case of arbitrary spin s, in which 
the size of the problem is reduced from 3(2s+l) x3(2s+l) 
down to (2s + 1) x (2s + 1), because the identity (|55j) 
projects the whole physical space onto the (0, 0)-charge 
sector. Of course, one has to pay the "price" for this re- 
duction of size: the appearance of BD~ 1 C term, in which 
information on the original larger matrix is compressed. 
The BD~ X G term is, as it should be, a square matrix of 
size (2s + 1) x (2s + l) (though matrices B and C are gen- 
erally not square). Each row and column of this matrix 
are characterized by a set of quantum numbers (s z ,s' z ). 
As a consequence of angular momentum selection rules, 
the matrix BD~ 1 C, or what we will call later X s or Y s j, 
(proportional to BD~ X C), becomes a tridiagonal matrix, 
i.e., only (s z , s z ), (s z , s z + 1) and (s z , s z — 1) elements are 
non-zero. 

The tridiagonal matrices, X s and Y Si j, have a very 
characteristic structure (see Appendix), which allows us 
to identify one of its eigenvectors by inspection, Eq. (|3~lj) 
or (|3"T|) . We will see that this eigenvector gives also a 
solution of Eqs. |26|) and |25|) . which turns out, acci- 
dentally, to be the right solution satisfying the condi- 
tion (|27p . The matrices X s and Y s j consist of Clebsch- 
Gordan coefficients at fourth order. The above charac- 
teristic feature of those matrices are actually due to some 
unconventional quartic relations among Clebsch-Gordan 
coefficients, which are summarized in the Appendix. 

The explicit matrix elements of Eq. ([5]) is given as 
four constituent block matrices A — D in Eqs. ((H])-©. 
These formulas contain Fermi distribution functions, and 
therefore, depend on the relative positions of dot levels 
with respect to the chemical potential of two leads. At 
zero temperature such formulas can be further specified 
in specific situations which apply when the temperature 
is small compared to the voltage bias: (i) two spin sectors 
in the bias window, (ii) only one spin sector in the bias 
window (see Fig. [T]). 



A. Two spin subsectors in the bias window 

We first consider the case of two spin subsectors in 
the bias window, fi^ > (e s +i/2, es-1/2) > MR- At zero 
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L 
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R 
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'7777777/) 



L 
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'77777777/ 



(iii) 



FIG. 1: The dot states and the Fermi window: (i) two spin 
subsectors, (ii) only one spin subsector (j — s — 1/2), in the 
bias window, (iii) similar to (ii), but here the "triplet" (j — 
s + l/2) subsector is in the bias window. The figure is drawn 
for s = 1/2. 



temperature, Eqs. ©-I© become 

A{s z ,s z ) = - E C ] { Jz ,s z ) = -2T L , 

j'=s±l/2 j*=-j',- ,3 

Bj(s z ,j z ) = T R \(s,s z ;l/2,a z \j,j z }\ 2 , 



C 3 (j z , S ' z ) = T L \(s,s' z ;l/2,a' z \j,j z )\ 2 , 
Dj(jz,jz) = - E B i( s ^iz) = -Ti 



The rate matrix ([5]) has the simple form: 



(29) 



B 

D + £I 



With the use of identity (j2"5|) , the eigenvalue problem 
defined by Eqs. (f2"rl)) and ([25]) reduces to finding A which 
satisfies, 

det [A -XI- B{D + (£ - X)I}- 1 Ce iri ] = 0. 

As mentioned before, the BD~ 1 C -like term is a (2s + l) x 
(2s + 1), tridiagonal matrix with finite matrix elements 
only at (s z ,s z ), (s z ,s z + 1) and (s z ,s z — 1). Recalling 
Eq. (USD, one finds, 



B{D + (£ — X)I}~ 1 Ce 



vr] 



T L T R e"' 



where X s is a (2s + 1) x (2s + 1) matrix, whose explicit 
matrix elements are, 

X s (s z ,s' z ) = Y, E l(s,^;i/2,^|i,i,)| 2 

j=s±l/2 jz=-j,~ ,3 

x \(s,s' z ;l/2,a' z \j,j z )\ 2 , (30) 
The reduced eigenvalue equation now reads, 

det [{-2T L X)I + ^fl^ X s ] = 0. 

As is explained in the Appendix, this matrix X s has the 
striking property: 



(31) 













1 




1 


X s 




= 2 






K 1 ) 




K 1 ) 



valid for arbitrary s. The vector vq — (1, 1, • • • , 1) T is 
thus an eigenvector of X s with an eigenvalue 2 for all s, 
but it is also an eigenvector of A — XI. This means that 
if A = A(£, ?y) is a solution of the quadratic equation, 



(32) 



(2r L + A)(r fl + A-p _ 2 

then satisfies, 

[A -XI- B{D + (f - A)J} -1 Ce*»]«o = 0, 

Consequently, this solution A is a solution of the origi- 
nal eigenvalue problem defined by Eqs. (f25|) and (|26|) . 
However, at this point it is still only "a" solution of the 
original eigenvalue problem. The remarkable feature of 
Eq. (|3"2")) is that in the limit of £, 77 — > 0, if one takes also 
the limit of vanishing counting fields: A — > 0, then on the 
left hand side, the numerator and the denominator cancel 
each other simply to give a numerical factor 2, which co- 
incides the right hand side of the equation. This implies 
that a solution of the quadratic equation (|3"2"|) is smoothly 
connected to the stationary state solution, which satisfies 
Eq. (|2T|) . Indeed, one of the solutions of Eq. 



2T L + T; 



iV(2r L -r fl + c) 2 + 8r L r^ (33) 



satisfies Eq. (|27|). This means that Eq. ([33]) is actu- 
ally the desired cumulant generating function yielding the 
FCS. Thus we are able to identify the suitable eigenvalue 
of the problem formulated as Eqs. (|2"rj)) and (|25[) simply 
by inspection. The obtained result (|33|) for general s is 
identical to the s — 1/2 case, because the corresponding 
eigenvalue of X s , i.e., Eq. (f3"Tj) is not dependent on s. 
Note that Eq. ([33)) is also equivalent to the the CGF cal- 
culated for spinful conduction electrons in the absence of 
local (molecular) spin. [l2j 



B. Only one spin subsector in the bias window 

Consider now the case of only one spin subsector, say, 
j = s — 1/2 in the bias window, i.e., e s+ i/2 > hl > 
e s-i/2 > fJ'R- This corresponds to the case (ii) in Fig. Q] 
In this case, Eqs. (|6])-((9j) reduce, at zero temperature, to: 



A(s z ,s z ) = -- 



2s 



2s + 1 

B s -i/ 2 (s z ,j z ) = T R \(s,s z ; 1/2, a z \j = s - l/2,j z )\ 2 , 

B s+ i/ 2 (s z ,j z ) = {T L +r R ) 

x \(s,s z ;l/2,a z \j = S + l/2, Jz )\ 2 , 

C.-1/aCj*,*;) - r L |(s, s' z ; 1/2, a' z \j = s - 1/2, j 

Cs+l/2(jz,s' z ) = 0, 

D s -i/ 2 (jz,j z ) = -Tr, 

D s +l/2(jz,j Z ) = -T L -T R , 



3/| 2 J 
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where u z and u' z are not independent variables, but they 
are determined automatically by the constraints: j z = 
s z + a z and j z — s' z + a' z , respectively. The modified rate 
matrix can be obtained as Eq. (j2"5)) . Again we use the 
identity, Eq. l[28]). Note that C s+1/2 = 0. This allows us 
to rewrite the BD~ X C term as, 



5.-1/2^-1/2^ 



/2 



r fl + a - e 



Y 



s,s— 1/2; 



where Y],j is a (2s + 1) x (2s + 1) square matrix, whose 
(s z , s' z )-elements are, 



YsA s ^ s 'z)= \(s,s z ;l/2,a z \j,j z )\ 2 

3m = —3t— ij 

x\( S> s' z ;l/2,a' z \j,j z )\ 2 . 
The eigenvalue equation reads, 



det 



2s 



2s + 1 



r L - a / 



r, 



A-e 



i — 1/2 



(34) 

= 0. 

(35) 

Another situation to be considered in parallel is the 
case of only j = s + 1/2 spin sector in the bias window, 
i.e., e s _i/2 > Ml > e s +i/2 > Mfl- This corresponds to the 
case (hi) in Fig. [T] In this case, Eqs. ©-(J!]) reduce, at 
zero temperature, to, 



A(s z , s z 



2s + 2 



2s+l 

B.-i /2 (a z ,j z ) = (r L + r fl .) 

x \{s,s z ;l/2,a z \j = s-l/2,j z )\ 2 , 
B s+ i/ 2 {s z ,j z ) = T R \(s,s z ;l/2 7 a z \j = s + l/2,j z )\ 2 , 

C s -l/2(jz,s' z ) = 0, 

C s+1/2 ( ]zi s' z ) = r L |(s,s' 2 ;l/2,a z |j = s + l/2,j z )| 2 , 
D s -i/2(jz,jz) = -T L -T R , 

D s +l/2(jzjz) = -Tr- 

The modified rate matrix can be constructed accordingly. 
Note that this time C s _i/2 = 0. After the use of identity 
}, the reduced eigenvalue equation reads, 



det 



2s + 2 
2s + I' 



Tr-X)I 



Y 



s,s+l/2 



= 0. 



(36) 



We now attempt to identify the maximal eigenvalue 
A(£, rj), following the same logic as the case of two spin 
subsectors in the bias window. As expected, the matrix 
Y s j has the following characterizing feature: 



Y 









1 


_ 2 j + 1 


1 




2s + 1 




w 




W 



(37) 



This is explained in the Appendix. Eqs. (|35|) . (136)) and 
(|37| suggests that A(£, 77) should satisfy a quadratic equa- 
tion analogous to Eq. (|32j) . 



(fSr 1 ^ + a) (r,, + a - 



2j + l 
2s + 1' 



(38) 



then A(£,77) satisfies also Eqs. (1331) and (|3"6"|) . Of course, 
at this point it is still only "a" solution of the original 
eigenvalue problem. However, Eq. (f3"2"| has again the 
following remarkable property: in the limit of vanishing 
counting fields, 77 — ^ 0, if one takes also the limit A — > 0, 
then on the left hand side the numerator cancels with the 
denominator to give simply a factor ^l+i j which is insen- 
sitive to the T's, and coincides with the right hand side 
of the equation. This a rather unexpected coincidence 
because the factor t^+T on the left hand side of the equa- 
tion comes from a quadratic relation (|A1[) , whereas the 
same ^l+i factor on the right hand side originates from 
a quartic relation (|A10|) . Indeed, one of the solutions of 
Eq. ([38]), 

A(Si , ) = _M^±£izi 




L 2j + 1 
2s + 1 



r L r R e^, 



(39) 



satisfies Eq. (|2"T1) . We are thus able to solve the problem 
simultaneously for the two cases, i.e., either of the spin 
subsectors (j = s=F 1/2) in the bias window, correspond- 
ing to the cases (ii) and (iii) in Fig. [TJ Contrary to the 
previous case, i.e., Eq. (|33|) f° r t w0 spin subsectors in 
the bias window, result which we obtain (|39[) depends on 
s. This is because the eigenvalue of Eq. (|3T|) is depen- 
dent on both s and j. The generating function A(£,?y) in 
Eq. (|39[) is also asymmetric with respect to the bare am- 
plitudes, and T R , whereas this asymmetry tends to 
disappear for large spin s. These characteristic features 
of the charge-current joint distribution function are fur- 
ther discussed in the next section. 



C. Remarks 

We have thus successfully calculated the FCS for trans- 
port through a molecular quantum dot magnet in the in- 
coherent tunneling regime. The results are obtained, in 
Eqs. ((33} and l[3"9")). together with Eq. $TZ§ in the form of 
a cumulant generating function (CGF) of charge-current 
joint correlation functions. Here, we apply the obtained 
CGF for deducing lowest-order current correlation func- 
tions, in order to check consistency with known results. 

In the case of two spin subsectors in the bias window, 
i.e., case (i) in Fig. [TJ the obtained CGF, i.e., Eq. (|3"3")1 
has no s-dependence. The lowest order cumulants are 
obtained by taking the derivatives of CGF with respect 
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to counting fields, and have no s-dependence. Thus the 
Fano factor is given for an arbitrary s by, 



Cq2 

C'oi 



ATI 



(2r z 



■r, 



(40) 



which coincides with the result for spinful conduction 
electrons and no local spin. [2lJ The skewness normalized 
by Coi reads, 

c 03 _ ier| - i6il + 24r|r| - at l t r + r% 



Cqi 



{2T L 



(41) 

Note that Eqs. (|4H|) and (|4"Tj) can be rewritten in terms 
of simple polynomials of a single parameter, x = {2Yl — 

r i? ) 2 /(4r L r iJ ) as, 



Co2 

Cqi 



1 



Cq3 

Cqi 



1 



Or + 2) 2 



One might expect that the (extended) Fano factors, 
C04/C01, C 05 /Coi, • ■ • could be obtained by developing 
this apparently systematic series. However, this naive 
expectation turns out to be an artifact at fourth order. 

In the case of j = s ±1/2 spin subsector in the bias 
window, one has to use Eq. (j3"3"j) . instead of Eq. ([55]) 
as the CGF. Notice that (|39|) is obtained by making the 



following replacement in Eq. ([55]) : 2Fl — » asTI^i,, e -£ 
the Fano factor becomes in this case, 



C02 
Cqi 



(any 



In the limit of large local spin s — » 00, this reproduces 
the result for spinless conduction electrons, 



C02 
Coi 



r 2 

1 R 



(42) 



The third order Fano factor (skewness normalized by 
Coi) is also obtained by applying the same replacement 
to Eq. gl]). 



V. DISCUSSIONS 

The exact analytic expressions for the CGF, Eqs. 
and (|39p . contain the full information of the non- 
equilibrium statistical properties of the current and of the 
charge. Though the CGF and the joint probability distri- 
bution contain the same information as we can see from 
their definitions (Sec. Ill), the former is rather a mathe- 
matical tool, whereas the latter provides us with a physi- 
cal picture. We therefore apply an inverse Fourier trans- 
formation to the obtained CGF, transforming it back to 
a joint probability distribution: 



P(N,I) 



(27T) 1 



drj / d£e' 



Q(£,ri) — iTN£-i tIt] 



The asymptotic form for r — ► 00 is obtained within 
the saddle point approximation, since Q is proportional 
to t, and thus the exponent is also proportional to 
t: In P(N, I) « Sl(C,V*) ~ iNr C ~ iIr V*, where f* 
and 77* are determined by the saddle point equations, 
NT = -id i n(C,V*) and I T = -id v Ct(£*,r)*) 12}. Com- 
bining these equations, we can produce a contour plot of 
In P(iV, /) as a parametric plot in terms of 77 and £. The 
probability distribution for only one of the observables, 
i.e., either charge or current, is obtained further by fixing 
77* — or £* =0, respectively. 

In this section, we try to visualize the joint probability 
distribution, and discuss how a molecular spin influences 
statistical properties of the SET molecular quantum dot. 
The joint distribution reveals nontrivial correlations: cor- 
relations among multiple current components in a multi- 
terminal chaotic cavity, [12j | or correlations between two 
different observables, such as current and charge, [l3| 
have been investigated. Here, we first consider a dot 
with symmetric tunnel junctions Tl — Tr, and analyze 
the joint distribution of charge and current under high 
and small bias voltages. We then consider the case of a 
large asymmetry in the tunnel coupling, e.g., Tl <C Tr. 
We will argue that much information on the molecular 
spin and the exchange coupling can be extracted from 
transport measurements. Finally, we will plot the lowest 
three cumulants, which are indeed being measured in the 
state-of-art experiments. 



A. Effects of local spin on the statistical properties 
for symmetric tunnel coupling Tl=Tr 

Let us consider a molecular quantum dot with sym- 
metric tunnel junctions = Tr. We first consider the 
case of a large bias voltage, in which the two spin sub- 
sectors are both in the bias window, i.e. the case (i) in 
Fig. [TJ In this case the CGF has been calculated ana- 
lytically in Eq. (|33p . Here, we attempt to uncover basic 
physical consequences, hidden in this result. In Figs. 
2(a) and 2(b), the probability distribution functions of, 
respectively, charge and current are plotted. They have 
been evaluated by applying the saddle-point approxima- 
tion explained above to the solution, Eq. (f3"3")l . In the 
two cases, the distribution function is not symmetric in 
regard to the peak position, indicating a clear deviation 
from the Gaussian distribution. Such asymmetric distri- 
butions were somewhat expected, since, e.g., in the case 
of charge distribution, the average value (N) is not lo- 
cated at 1/2 but at 2/3. 

In Fig. 2(c) the two distributions are shown simulta- 
neously as a joint distribution in the form of a contour 
plot. The distribution possesses a single peak at 1= (I) 
and N = (N) . For a given value of current I, the peak 
position of N converges to N = (N) = 2/3 as the current 
becomes larger (I^> (I)). Note that the probability dis- 
tributions depend neither on the size s of the molecular 
spin nor the exchange coupling g. The asymmetry in the 
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FIG. 2: The probability distributions of charge (a) and cur- 
rent (b) for the symmetric coupling Tl — Tr. The bias voltage 
is large enough and both of the two spin subsectors are in the 
bias window [Fig. [1] (i)]. (c) The contour plot of the non- 
equilibrium charge-current joint distribution \nP(N, I). The 
horizontal (vertical) axis corresponds to current / (charge N). 
The contour interval is (2T L +Tr)t /20. 



charge distribution comes from the coefficient 2 in front 
of [Eq. l[55])]. reflecting the number of degrees of free- 
dom of a transmitted spin. The truth is that Eq. (|33D is 
not different from the CGF in the absence of a molecular 
spin, i.e., the coefficient 2 appears also in the latter case. 
[r2j | We thus conclude that at a high enough bias voltage, 
all the statistical properties will reproduce the results in 
the absence of a molecular spin. 

We then turn our attention to the small bias voltage 
regime, where only one spin subsector is in the bias win- 
dow, the cases (ii) or (iii) in Fig. [TJ Contour plots of the 
joint probability distribution obtained by using Eq. (|39|) 
are plotted in Fig. [3] Figs. 3(a) and 3(b) correspond 
to a small molecular spin s = 1/2. In both ferromag- 
netic (g < 0) [panel (a)], and antiferromagnetic (g > 0) 
[panel (b)] cases, the charge distribution is not symmet- 
ric, i.e., typically the peak position is away from N = 1/2. 
The nature of this asymmetry is similar to the previous 
case, i.e., the case of large bias voltage, Fig. [2] (c). The 
only difference is that here the peak position is located 
at (N) = ^T L /(^T L +T R ), which is 3/5 for panel (a) 
and 1/3 for panel (b). 

On the other hand, the distribution becomes symmet- 
ric w.r.t. a horizontal line N = 1/2 in the limit of a 
large spin s [panel (c)]. It might be a little surprise that 
this limit reproduces the joint probability distribution of 



a spinless fermion [l3[. In physical terms this can be 
interpreted as follows: a large molecular spin acts as a 
spin bath for the electron spin transmitted through the 
quantum dot. The n — charge sector is (2s + l)-fold 
degenerate, whereas the degeneracy of n — 1 charge sec- 
tor is either 2s or 2s + 2. For a large molecular spin s, 
the conduction electron spin does not feel the difference 
between the two different total angular momentum sub- 
sectors. The above result is counterintuitive in the sense 
that naively we expect that a large molecular spin may 
act as a large effective Zeeman field for transmitted spins. 




FIG. 3: Contour plot of the joint probability distribution 
\nP(N, I) in the case of only one spin subsector in the bias 
window [Figs. [T] (ii) and (iii)]. The coupling to reservoirs 
is assumed to be symmetric, Tl = Tr. Panels (a) and (b) 
correspond to a small local spin s — 1/2, whereas panel (c) 
corresponds to a large spin s — > oo. The exchange coupling 
between the local and an itinerant spin is ferromagnetic (g < 0) 
in panel (a), whereas antiferromagnetic (g > 0) in panel (b). 
The contour interval is (ffnTi +Tr)/20. 
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B. Asymmetry in the tunnel coupling Yl^Yr 

Here we discuss how an asymmetry in the tunnel cou- 
pling may change the statistical properties of a SET 
molecular quantum dot. In reality, the tunneling cou- 
pling to reservoirs is unlikely to be symmetric, because 
one cannot control the coupling strength of ligands to 
metallic leads @. In the extremely asymmetric limit 
;§> Tl, the CGF of current reduces to a Poissonian 
form: 

fi«Tzr L (e" ) -l). (43) 

The factor z is a positive constant, which is smaller than 
2: z = 2 for the case (i) and z = (2j + l)/(2s + l) with 
j = s=Fl, respectively, for the cases (ii) and (iii) in Fig. [J 
When we apply a negative bias voltage, the CGFs are 
obtained from Eqs. (l33|) and ([39]) by replacing Yl^Yr 
and 77 — — 77 as, QktYl (e~ tTI — 1). The absolute value 
of the n-th cumulant is |Co ra | ~tzYl for a positive bias 
voltage, whereas it becomes \Co n \ ~t~Yl for a negative 
bias voltage. Therefore, from the ratio |Co n /Conl> we 
can estimate the factor z, and consequently, obtain some 
information on the molecular spin and the exchange cou- 
pling. Such a method was previously proposed [22| | for a 
multi-orbital quantum dot. 

Figs. 4(a) and 4(b) are the joint probability distri- 
bution for Tl/Tji = 0.1. A strong asymmetry around 
N={N) = zT L /(zT L + T R ) is observed both for the an- 
tiferromagnetic (a) and the ferromagnetic (b) cases. A 
longer tail in the panel (b) implies that the correlations 
among tunneling processes are weaker for antiferromag- 
netic coupling. The measure of such correlations, the 
Fano factor C02/C01, as a function of the asymmetry ra- 
tio Tl/Yr is shown in Fig. [4] (c). We can observe that 
even for a small asymmetry ratio Y l/Y r = Q.1, the Fano 
factor is still smaller than 1 , and a fair amount of differ- 
ence remains between the two cases. 

For the charge distribution, around / « 0, we observe 
equidistant contours in panels (a) and (b) , which implies 
that the charge is exponentially distributed. Actually 
the charge distribution for both cases roughly follows 
In P(iV) oc — tTrN [Figured] (d)]. A similar exponen- 
tial distribution appears as an equilibrium distribution 
of charge when the dot level is far away from either of 
the two lead chemical potentials [l3j]. It is shown that 
the exponent depends simply on the decay rate of the 
excited state. In the present case, since the stationary 
state is fairly approximated by the empty state, this re- 
duces to the outgoing tunneling rate Yr of an electron in 
the occupied state to the right reservoir. 



C. Relation to experiments 

Up to now, we have shown contour plots of the proba- 
bility distribution of charge and current. For a semicon- 
ducting quantum dot, it has become possible to measure 




0.5 1 

n/r R 




N 



FIG. 4: Contour plots of InP for a small local spin s = 1/2 
for the asymmetric coupling Yr = 10IT. The ferromagnetic 
(g < 0) and antiferromagnetic (g > 0) coupling cases are shown 
in panels (a) and (b), respectively. The contour interval is 
(r_L+r_R)/20. (c) Fano factor as a function of the asymmetry 
factor Yl/Yr and (d) Probability distributions of charge for 
antiferromagnetic and ferromagnetic couplings. 

the statistical probability distribution of current in the 
incoherent tunneling regime [j"?], EH • Such measurements 
are based on the real-time counting of single electron 
jumps using an on-chip quantum point contact charge de- 
tectors. At the present stage, it may be rather challeng- 
ing to apply the same technique for a single-molecular 
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device with the same precision, but it might be still pos- 
sible to measure the lowest three cumulants, since the 
skewness for a tunnel junction can be measured by the 
state-of-art experiments. [13 . [l5| 

We have shown explicitly the lowest order cumulants 
in Sec. IV C. For a large bias voltage, they are given, e.g., 
in Eqs. (PD|) and PT]) . For a small bias voltage, the (gen- 
eralized) Fano factors were obtained by making the re- 
placement, 2T L -> zT L = §j^j-r L , in Eqs. (0OD and (gTJ). 
Fig. 5 shows the bias voltage dependence of the lowest 
three cumulants for a large molecular spin s = 10, which is 
approximately the size of the magnetic moment of Mni2. 
@ In this case, as the factor z = = 20/21 (for 

j = s — 1/2) is close to 1, the molecular spin is expected 
to act as a spin bath, reproducing characteristic features 
of the transport of spinless fermions. It is assumed that 
the tunnel coupling is asymmetric, = 5Tl. We also 
focus on the antiferromagnetic (g > 0) case, in which the 
"triplet" level (j = s + 1/2), e s+1 / 2 = edot + ffs/2 is above 
the "singlet" level (j = s-1/2), e s _ 1/2 = e dot -g(s+l)/2. 
In equilibrium, the singlet subsector is tuned to be at 
the lead chemical potential level, /il = £ s -\/2 = = 0. 
We then apply either a positive [panels (a-1) and (b- 

1) ] or a negative [panel (a-2) and (b-2)] bias voltage. 
The singlet level always stays in the bias window, i.e., 
either [i^ > e s _i/ 2 > \ir [panels (a-1) and (b-1)] or 
Hl < (-s-1/2 < [panels (a-2) and (b-2)] is realized. 
Note that in Fig. 5 the origin of energy is taken, for sim- 
plicity, to be at the singlet level, i.e., e s ~i/2 = is always 
satisfied. 

The expected cumulant-voltage characteristics for a 
positive and for a negative bias voltage are depicted 
in panels (a-1) and (a-2) of Fig. 5. At a small bias 
voltage, \hl,r\ < e s +i/2 - e«-i/2 = 9( s + 1 / 2 ), the ab- 
solute value of various cumulants are insensitive to the 
sign of bias voltage [This case realizes the situation of 
Fig. [2] (c)]. On the other hand, when the bias voltage is 
large, \hl,r\ >5( s + 1/2), each cumulant takes different 
values depending on whether the bias is either positive 
\fiz,\ > g(s+l/2), or negative \hr\ > g(s + 1/2) [This case 
realizes the situation of Fig. [3] (c)] . Such a feature could 
be a smoking gun of the transport of spinless fermions. 
Panels (b-1) and (b-2) show Fano factors corresponding, 
respectively, to the panels (a-1) and (a-2). At a large 
bias voltage \hl,r\ >.9( s + 1/2), the suppression of Fano 
factors is stronger for a positive [panel (b-1)] than for a 
negative [panel (b-2)] bias voltage. This can be under- 
stood as follows. If one compares the large bias regime of 
panels (a-1) and (a-2), one notices that the absolute value 
of current, i.e., Coi is smaller in this regime in panel (a- 

2) . This implies that tunneling events are less correlated 
when the bias is large and negative, leading to smaller 
Fano factors in panel (b-2). This is rather a trivial fact, 
but subjective to a direct experimental check. 

Finally, it should be noted again that the magnetic 
anisotropy of the molecule was neglected throughout the 
present paper. The extension of our theory to account for 
a large anisotropy comparable to the exchange coupling 
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FIG. 5: The absolute value of the lowest three cumulants, 
Coi, C'o2 and C03, for a molecular spin s — 10 and an asym- 
metric tunnel coupling (r.R = 5r£,). We also assumed that the 
exchange interaction is antiferromagnetic (g>0). The origin 
of energy is taken to be at the "singlet" level, e s _i/2 = 0. 
Panels (a-1) and (a-2) corresponds, respectively, to a positive 
and to a negative bias voltage. In panels (b-1) and (b-2), the 
corresponding Fano factors C02/C01 and C03/C01 are plotted. 



is beyond our scope here. 



VI. CONCLUSIONS 

We have thus studied the full counting statistics (FCS) 
for a molecular quantum dot magnet, a quantum dot with 
an intrinsic molecular spin (s) degrees of freedom, which 
is coupled to the conduction electron spin (er) via an ex- 
change interaction (s • a). Such a molecular quantum 
dot magnet is now experimentally available, and carries 
generally a large molecular spin s ~ 10. [6| We applied 
the master equation approach to FCS 12] to our molecu- 
lar quantum dot problem, by considering the incoherent 
tunneling regime. In this approach the cumulant gen- 
erating function (CGF) for FCS is obtained by solving 
an eigenvalue problem associated with a modified mas- 
ter equation with counting fields. We also assumed that 
the molecular quantum dot was in the strong Coulomb 
blockade regime (U — > 00) so that the dot state has 
only two charge sectors n = 0,1. A standard algebraic 
identity, Eq. (|28p. allowed us to reduce the size of the 
eigenvalue problem from 3(2 s + 1) x 3(2 s + 1) down to 
(2s+l) x (2s+l) by projecting the relevant physical space 
onto the charge (0, 0)-sector. This was, of course, simply 
a rewriting of the original problem. In the new repre- 
sentation, Clebsch-Gordan coefficients appear at fourth 
order. We established a few identities among such fourth 
order coefficients, which allowed us to identify one solu- 
tion of the eigenvalue problem, which turned out to be 
the correct eigenvalue, representing the CGF for FCS. 
Thus we obtained an analytic expression for the CGF as 
a function of the local spin s for different configurations 
of the dot states and of the two leads in regard to their 
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relative energy levels (Fig. [J). 

Based on the obtained analytic expressions, we also 
developed numerical analysis, e.g., contour plots of the 
charge-current joint distribution function. For a small 
local spin, e.g., s = 1/2, the obtained contour plots show 
clearly an asymmetry in the distribution of the charge, 
i.e., between n = and n = 1. We then demonstrated 
that this asymmetry grows differently for different cases 
in Fig. [1] mentioned above. In particular, this asymme- 
try tends to disappear as the local spin s becomes large, 
in cases (ii) and (iii) of Fig. [TJ indicating that the local 
spin plays the role of a spin bath for the conduction elec- 
tron. Of course, as emphasized in Sec. 5, it is rather 
counterintuitive to recover such characteristic transport 
properties of a spinless fermion in this large s limit. We 
further demonstrated that such characteristic features of 
the transport of a spinless fermion, which we believe to 
have a chance to be realized in reality, e.g., in a Mni2 
molecular quantum dot, are subjective to a direct ex- 
perimental check through the bias voltage dependence of 
lowest-order cumulants. 

It is interesting to extend our analysis to the case of 
leads of different nature, such as ferromagnetic leads or 
superconducting leads. The range of validity of this work 
is limited to the incoherent tunneling regime, whereas one 
might naturally wonder what happens when the coupling 
between the dot and the leads become stronger. Kondo 
typ e physics in such a regime is particularly interesting, 
|23| because of the coexistence of an intrinsic magnetic 
impurity (molecular spin) and the quantum dot, the lat- 
ter known to show enhancement of conductance due to 
Kondo mechanism at low temperatures. Certainly one 
has to go beyond the master equation approach to FCS 
in order to treat FCS in such a regime. 
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Appendices 



APPENDIX A: USEFUL RELATIONS AMONG 
CLEBSCH-GORDAN COEFFICIENTS 

For the calculation of FCS for a model with an arbi- 
trary spin s we have used extensively the Clebsch-Gordan 
algebra. We, therefore, summarize in the following useful 



relations among Clebsch-Gordan coefficients. We start 
with conventional quadratic relations which are some- 
times discussed in standard quantum mechanics text- 
books, and then proceed to unconventional quartic rela- 
tions which we discovered as a byproduct of this research. 



1. Quadratic relations 

Let us summarize the quadratic identities on the 
Clebsch-Gordan coefficients, (s,s z ;l/2,a z \j,j z ), for a 
given s and a = 1/2. Angular momentum selection rule 
allows only two values of j — s ± 1/2. Because of the 
constraint j z — s z +a z , only two of the three parameters, 
s z , a z and j z are independent. Lets us first look at Ta- 
ble I. The two rows correspond to different total angular 
momentum j. Either of the two columns specified by o~ z 
correspond also to the same j z = s z +a z . Then summing 
up the two coefficients squared given in either of the two 
rows, one finds, 

\( S , Sz ;l/2,a z \j, Jz )\ 2 = 

o»=±l/2 

J2 \(s,s z ;l/2,a z \j, Jz )\ 2 = g±i. (Al) 

3»=-jV" i3 

The identity (|A1[) applies to a given set of s, j and s z . If 
we release the quantum number j in Eq. IjAip , and take 
the summation over two spin sectors, one finds, 

J2 E \(s,s z ;l/2,a z \j,j z )\ 2 

j=*±l/2 U=-3,— ,3 

2s 2s + 2 „ 

= 1 = 2. A2 

2s + 1 2s + 1 v ' 

which is independent of the spin s. The identi- 
ties (|Aip , (|A2[) are used to obtain the expressions for 
A(s z ,s z ). 





a z = -1/2 


o z = 1/2 


j = 8- 1/2 


V 2s + l 


/ s ~ s z 
V 2^+1 


3 = * + 1/2 


/«-»,+! 
V 2s+l 


/s + s.+l 
V 2s+l 



TABLE I: Clebsch-Gordan coefficients, {s,s z ;l/2, a z \j, jz = 
Sz+Cz) for given s and s z . Note the difference between Tables 
I and II. They are equivalent but represented in different ways 
for different use. The four elements on the Table correspond 
one by one, i.e., they are same in number in different Tables, 
but written in terms of different parameters (either j z or s z ). 

On the other hand, keeping s and j fixed, one can also 
consider another situation where j z is fixed (instead of 
s z ). This corresponds to Table II. If one sums up the 
two coefficients squared given in either of the two rows 
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o z = -1/2 


(7 2 = 1/2 


j = 8- 1/2 


V 2s + l 


/«-J*+l/2 
V 2s + l 


3 = s + 1/2 


/s-ia + 1/2 
V 2s + l 


/s+jz + 1/2 
V 2s + l 



TABLE II: Clebsch-Gordan coefficients, (s, s z = j z — 
a z ;l/2,a z \j,j z ) for given s and j z . Note that for s and j z 
given, only four elements shown in the Table remain finite. 



in Table II (as we did for Table I), one finds this time, 
53 \(s, Sz ;l/2,a z \j,j z )\ 2 = 

CT=±l/2 

53 |<s,s z ;l/2,a z |j,j 2 }| 2 = 1. (A3) 



s, — — s, ■■■ ,s 



If we focus on a particular term in the above summation, 
for which s z is given, then a z is determined automatically 
by the constraint j z = s z + a z . Therefore, the summa- 
tion over <7 2 and s z can be replaced by one another. The 
identity (|A3|) is relevant for finding the diagonal elements 
of the matrix D defined by Eq. ©. It turns out that 
Dj (j z , j z ) does not depend on j z . The identity (|A3|) ap- 
plies to a given set of s, j and j z . 

The two identities (|A1IA3[) are indeed consistent, 

53 53 \(s,s z ;l/2,a z \j,j z )\ 2 = 

2.7 + 1 



E 



2 J + 1, 



2s + 1 

53 53 \(s,s z ;l/2,a z \j,j z }\ 2 = 

53 l = 2j + l. (A4) 



Jz=~3,--- ,3 Sz=—8,— ,s 



3*=-lt— ,3 



Furthermore, releasing again the quantum number j and 
taking the summation over two spin sectors, one can ver- 

ify, 

53 53 53 |(s,s 2 ;l/2,a 2 |j, j2 )| 2 

s z =-s t - ,sj=a±l/2j*=-j,- ,j 

= 53 2 = 2(2s + l), 



S z — — S,'" ,s 



E E E \(s,s z ;l/2,a z \j,j z }\ 2 

j=s±l/2js=-j,— ,3 »«=-»,— ,s 

= 53 (2j + l) = 2(2s + l). 

We used Eq. (|A2|) for obtaining the second expression, 
whereas for the third one, we used Eq. (|A4[) . but, of 
course, whichever path one chooses, one always obtain 
at the end, the same 2s + 1. 



2. The quartic relation 

As a by-product of the calculation of FCS, we dis- 
covered some unconventional quartic relations among 
Clebsch-Gordan coefficients. In the evaluation of a deter- 
minant derived from the modified master equation (|23p 
we encountered a (2s+l) x (2s+l) square matrix, consist- 
ing of forth-order Clebsch-Gordan coefficients, defined as 
Eq. (|3"0|) , whose (s z , s 2 )-elements one can rewrite here 
as, 

x s (s z ,s' z )= 53 E 

j=a±l/2j„=-j,- ,3 

\(s,s z ;l/2,a z \j,j z )\ 2 \(s,s' z ;l/2,a' z \j,j z )\ 2 . 

X s is a tridiagonal matrix, i.e., only (s z , s z ), (s z , s z + 1) 
and (s 2 ,s 2 — 1) elements are non-zero. They are given 
explicitly as, 

X s (s z ,s z ) 

X s (s z ,s z + 1) 

*.(..,..-!) - 2(s+ S + ;; +1) <-» 

This matrix X s has a very characteristic property, Eq. 
@D, i.e., 



2{s 2 


+ (s + l) 2 - 


^2s 2 } 




(2s +1) 2 




2(a- 


s«)(s + « 2 


+ 1) 




(2s + 1)2 


5 


2(s + 


s 2 )(s - s 2 


+ 1) 


(2s + 1) 2 



V 1 / V 1 / 

This can be checked explicitly for s 
Indeed, the matrix X s reads, 



1/2,1,3/2,- 



X 



1/2 



/ 3/2 1/2 
I 1/2 3/2 



X- 



3/2 



/ 14/9 4/9 

,Xi= 4/9 10/9 4/9 

y 4/9 14/9 

f 13/8 3/8 \ 
3/8 9/8 1/2 
1/2 9/8 3/8 
V 3/8 13/8 j 



(A6) 



Proving the general result (|A6|I for arbitrary s is equiva- 
lent to establishing the following identity: 

E E E 

j=s±l/2 «i=— s,~ ,s ]*=—],••• J 

K S ,s z ;l/2,a z |j ! i,)| 2 |( S ,4;l/2,^|j ! i,)| 2 = 2. (A7) 

When s 2 ^ ±s, one can verify this by substituting into 
it explicit matrix elements given in Eq. (|A5[) . i.e., 



X s (s 2 , s, - 1) + A s (s 2 , s 2 ) + X s (s z , s z + 1) = 2. (A8) 
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Then one can further verify that Eq. (|A8|) still holds 
when the edges of the X s matrix are encountered, at s z = 
±s. 

The same kind of structure as Eqs. (|A6|) and (|A7|) 
exists actually for each spin sector j. Y s j introduced in 
Eq. ()34|) is such a matrix, whose (s z , s' z )-elements one 
can rewrite here as, 

Y s ,j(s z ,s' z ) = ^2 \(s,s z ; 1/2, a z \j, j z )\ 2 
ix=-j>— J 

x\(s,s' z ;l/2,a' z \j,j z )\ 2 . 

Y sJ is related X s as, X s = J2j= s ±i/2 Y,j = Y s,s-i/2 + 
Y s , s +i/2- Again Y s _j is a tridiagonal matrix, with fi- 
nite matrix elements only at (s z ,s z ), (s z ,s z + 1) and 
(s z , s z — 1). As expected, the matrix Y s j has the follow- 
ing characterizing feature, Eq. (|57|) . i.e., 



(l\ 



2j + l 
2s + 1 



V 1 / 



(A9) 



which is also equivalent to the following identity: 
]T Y, \{s,s z -A/2,a z \j,j z )\ 2 

2j + 1 



x\(s,s' z ;l/2,a' z \j,j z )\ 2 



2s + 1 



(A10) 



Now in order to prove Eq. (|A9|) or (| A10[l for general s 
let us consider the following two cases separately: 

1. j = s — 1/2 spin sector. One can first directly verify 
Eq. (|37[) by writing down explicitly the matrices 
n,.-i/2 for 5 = 1/2,1,3/2, 



Yi 



/2,0 



1/2 



3/2,1 



1/4 1/4 \ 
1/4 1/4 J ' 

4/9 2/9 
2/9 2/9 2/9 
2/9 4/9 

/ 9/16 3/16 \ 

3/16 5/16 1/4 

1/4 5/16 3/16 

\ 3/16 9/16 / 



The finite tridiagonal matrix elements are in this 
case, 



Y. 



i-l/2 



(Sz,S z ) 



2(s 2 



Y 



s.s — l 



/2(S Z ,S Z + 1) 



(2s + l) 2 ' 
(s - s z )(s + s z + 1) 
(2s + 1) 2 



v I -n (s + Sz)(s-s z + l) 
Y s , s ^ /2 (sz,s z -l) = ( 2s + i)2 -( AU ) 



2. j = s + 1/2 spin sector. One verifies Eq. (137)) by 
writing down explicitly the matrices Y s s +i/2 f° r 
5 = 1/2,1,3/2,-..: 



^1/2,1 
Y,S/2 



Yp 



3/2 : 2 



5/4 1/4 \ 
1/4 5/4 J ' 

10/9 2/9 
2/9 8/9 2/9 
2/9 10/9 

/ 17/16 3/16 \ 

3/16 13/16 1/4 

1/4 13/16 3/16 

V 3/16 17/16 J 



Notice also that X s = Y StS _i/2 + F s s+1 / 2 . The 
finite tridiagonal matrix elements are, 



Y s<a +l/2\ s Zi Sz) 

Y SiS+1 / 2 (s z , s z + 1) 
Y StS+ i/2(s z , s z - 1) 



2{(s + l) 2 + ^} 
(2s + 1)2 ' 

(a - s z )(s + s z + 1) 

(2s + 1)2 
(s + s z )(s-s z + l) 

(2s + 1)2 



(A12) 



When s z ^ ±s, one can verify Eq. (|A9[) . using either Eq. 
(|A11|) or (|A12|) . i.e., one finds for a given s z , 



Y 



S,J \°Z 1 



l) + Y a<j (s z ,s M ) 



+Y Stj (s z ,s z + 1) 



2j + l 
2s + 1' 



(A13) 



Then one can further verify that Eq. (|A13|) still holds 
when the edges of the X s matrix are encountered, at s z = 
±s. This establishes Eq. (|A"9j) . or equivalently, (|A10|> . 



Note that the quartic summation, Eq. (|A10[) . or the 
eigenvalue of Eq. (|A9|) happens to be identical to a 
quadratic summation, (|A1[) . This coincidence is a crucial 
ingredient for the FCS solvability of our model. If one 
takes the summation over two spin sectors j = s ± 1/2, 
Eq. (|A10|) recovers Eq. ([AT]) . 
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